% Influence, by size decile
bins=3;
lambda = zeros(bins,1);
lambda_h = zeros(bins,1);
line = 1;
for i = [90 9 1]
    fn = strcat('MAT/rho',num2str(rhoT),'eta',num2str(etaT),...
        'lambda',num2str(lambdaT),'_psi',num2str(varphi),...
        '_sigma',num2str(sigmaT),'/Step_6_output_alpha',...
        num2str(alphaT),'_rho',num2str(rhoT),'_DOMidioprod_approx_Percentile',...
        num2str(i));
    load(fn,'e_Y_FRf_doubledef');
    lambda(line) = e_Y_FRf_doubledef;
    
    fn = strcat('MAT/rho',num2str(rhoT),'eta',num2str(etaT),...
        'lambda',num2str(lambdaT),'_psi',num2str(varphi),...
        '_sigma',num2str(sigmaT),'/Step_6_output_alpha',...
        num2str(alphaT),'_rho',num2str(rhoT),'_DOMidioprod_approx_Hom_Percentile',...
        num2str(i));
    load(fn,'e_Y_FRf_doubledef');
    lambda_h(line) = e_Y_FRf_doubledef;
    
    line = line+1;
end

% Exposure, at firm-level
eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',sigmaT,'/Step_4_setup_alpha',num2str(alphaT),...
    '_rho',rhoT,'_approx.mat')]);
eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',sigmaT,'/Step_5_output_alpha',num2str(alphaT),...
    '_rho',rhoT,'_WORLDprod_p10','_approx.mat')]);
eval([strcat('load MAT/WIOD.mat')]);
eval([strcat('load MAT/Step_2_firm_data_alpha4.mat')]);

a_hat = ones(J*N,1);
a_hat = a_hat/1.1;
a_hat(start(15):start(15+1)-1) = 1;
g_A = a_hat-1;
I = eye(J*N);
K = transpose(alpha_WIOD_nj);
K = repmat(K,[1,J*N]);
B = transpose(gamma);
A = I-B.*K;
g_P = inv(A)*g_A;
P_hat = g_P +1;
pi_M_f_foreign = pi_M_f;
%pi_M_f_foreign(:,start(france):start(france+1)-1) = 0;
pi_M_f_foreign_g_P = pi_M_f_foreign.* repmat(transpose(g_P),[FI_J,1]);
sum_pi_M_f_foreign_g_P = sum(pi_M_f_foreign_g_P,2);
one_minus_alpha = 1-pi_l_f;
one_minus_alpha_sum_gamma_foreign_g_P = one_minus_alpha .* sum_pi_M_f_foreign_g_P;
exposure = one_minus_alpha_sum_gamma_foreign_g_P;
eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',sigmaT,'/Step_4_setup_alpha',num2str(alphaT),...
    '_rho',rhoT,'_approx.mat')]);
a_hat = ones(J*N,1);
a_hat = a_hat/1.1;
a_hat(start(15):start(15+1)-1) = 1;
g_A = a_hat-1;
I = eye(J*N);
K = transpose(alpha_WIOD_nj);
K = repmat(K,[1,J*N]);
B = transpose(gamma);
A = I-B.*K;
g_P = inv(A)*g_A;
P_hat = g_P +1;
X_nkif = firmsecD_sorted*X_mnj(start(france):start(france+1)-1,:);
pi_X_nkif=pi_mnjf.*X_nkif;
pi_X_nif=sum(pi_X_nkif,2);
one_minus_pi_l_f=1-pi_l_f;
pi_X_laborshare_nif=pi_X_nif.*one_minus_pi_l_f;
clear pi_X_nkif  X_nkif
pi_X_laborshare_nif_repmat = repmat(pi_X_laborshare_nif, [1 J*N]);     
gamma_pi_X_laborshare_nif = pi_X_laborshare_nif_repmat.*pi_M_f;
% Sum over firms. Should give a J*N matrix, which is numerator of
% homogeneous (1-alpha)*gamma
numerator = firmsecD_sorted'*gamma_pi_X_laborshare_nif;
% Sum exports over destination country per sector for france for numerator
% 32*1 matrix, which must be expanded per destination market
X_fra = X_mnj(start(france):start(france+1)-1,:);
den = sum(X_fra,2);   
denominator = repmat(den,[1 J*N]);
clear pi_X_laborshare_nif_repmat gamma_pi_X_laborshare_nif X_fra 
% Create homogeneous (1-alpha)gamma
one_minus_alp_gam_hom = numerator./denominator;
clear numerator denominator;
% Create homogeneous alpha
pi_X_laborshare_nif=pi_X_nif.*pi_l_f;
numerator = firmsecD_sorted'*pi_X_laborshare_nif;
pi_l_hom = numerator./den;
clear numerator den pi_X_laborshare_nif;
% Create homogeneous gamma given that alpha_hom = alpha_baseline at
% country-sector level
pi_l_hom_check = pi_l(start(france):start(france+1)-1,:);
corr(pi_l_hom,pi_l_hom_check) 
pi_M_hom = one_minus_alp_gam_hom./(1-repmat(pi_l_hom,[1 N*J]));
% Now fill in firm levels given homogeneous values and sectors firms in
pi_l_f_hom = firmsecD_sorted*pi_l_hom;
pi_M_f_hom = firmsecD_sorted*pi_M_hom;
clear denominator numerator one_minus_alp_gam_hom pi_l_hom pi_M_hom;
pi_M_f_foreign = pi_M_f_hom;
%pi_M_f_foreign(:,start(france):start(france+1)-1) = 0;
pi_M_f_foreign_g_P = pi_M_f_foreign.* repmat(transpose(g_P),[FI_J,1]);
sum_pi_M_f_foreign_g_P = sum(pi_M_f_foreign_g_P,2);
one_minus_alpha = 1-pi_l_f_hom;
one_minus_alpha_sum_gamma_foreign_g_P = one_minus_alpha .* sum_pi_M_f_foreign_g_P;
exposure_h = one_minus_alpha_sum_gamma_foreign_g_P;

% Aggregate across firms within a size quantile
eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',sigmaT,'/Step_5_output_alpha',num2str(alphaT),...
    '_rho',rhoT,'_WORLDprod_p10','_approx.mat')]);

ps=1;
PS=[];
for j=1:J
    perc_sec = zeros(size(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)));
    p = prctile(full(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)),99);
    perc_sec(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)>=p & perc_sec ==0) = 1;
    p = prctile(full(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)),90);
    perc_sec(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)>=p & perc_sec ==0) = 9;
    perc_sec(perc_sec ==0) = 90;
    PS = [PS;perc_sec];
end

exposure_bin = zeros(1,bins);
exposure_h_bin = zeros(1,bins);
line = 1;
for i = [90 9 1]
   exposure_bin(1,line)= mean(exposure(PS==i));   
   exposure_h_bin(1,line)= mean(exposure_h(PS==i));   
   line=line+1;
end

ratio_lambda = lambda./lambda_h./sum(lambda).*sum(lambda_h);
ratio_exposure = exposure_bin ./exposure_h_bin;

figure_8_line([1 2 3],ratio_lambda,ratio_exposure)
fig8 = strcat('Figures/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Fig8.eps');
saveas(gcf,fig8,'epsc2')


